unit Projection;
{************************************************************************}
{*     Projection.p                                                                                                                                          *}
{*     by Michael Castle (Pascal) and Janice Keller (assembly code)                                                          *}
{*     University of Michigan Mental Health Research Institute (MHRI)                                                     *}
{*     e-mail address: mike.castle@med.umich.edu                                                                                       *}
{* ** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *   * * * * * * * * * }

interface

	uses
		QuickDraw, PrintTraps, Palettes, Globals, Utilities, File2, File1, Graphics, Camera, Filters, Stacks;


	procedure DoProjection;
	function ShowProjectDialogBox: boolean;
	procedure Project;


implementation

	const
		bigpowerof2 = 8192;                               {used for integer trigonometric arithmetic}

	type
		MHRIRealArray = array[1..MaxPics] of real;
		RealPoint = record
				x: real;
				y: real;
			end;     {record}
		IntPtr = ^integer;
		LongPtr = ^LongInt;

	var
		SliceInterval: real;                                 {distance, in pixels, between original slices}
		BoundRect: rect;                                      {boundary rectangle for a rectangular selection}
		cancelled: boolean;
		lower, upper: integer;
		ProjSize: LongInt;



{******************************************************************************}
{*     This procedure frees memory allocated for buffers used in projection calculations.                                       *}
{******************************************************************************}
	procedure DisposeProjectionPtrs (projptr, opaptr, brightcueptr: ptr; zbufptr, countbufptr, cuezbufptr: IntPtr; sumbufptr: LongPtr);
	begin
		if zbufptr <> nil then
			DisposPtr(Ptr(zbufptr));
		if opaptr <> nil then
			DisposPtr(opaptr);
		if sumbufptr <> nil then
			DisposPtr(Ptr(sumbufptr));
		if countbufptr <> nil then
			DisposPtr(Ptr(countbufptr));
		if cuezbufptr <> nil then
			DisposPtr(Ptr(cuezbufptr));
		if brightcueptr <> nil then
			DisposPtr(brightcueptr);
		if projptr <> nil then
			DisposPtr(projptr);
	end;


{******************************************************************************}
{*     This procedure projects each pixel of a volume (stack of slices) onto a plane as the volume rotates about    *}
{*  the x-axis.  Integer arithmetic, precomputation of values, and iterative addition rather than multiplication  *}
{*  inside a loop are used extensively to make the code run efficiently.  Projection parameters stored in global   *}
{*  variables determine how the projection will be performed.  This procedure returns various buffers which   *}
{*  are actually used by DoProject to find the final projected image for the volume of slices at the current angle.*}
{******************************************************************************}
	procedure DoOneProjectionX (nSlices, ycenter, zcenter: integer; projwidth, projheight: LongInt; costheta, sintheta: longint; projptr, opaptr, brightcueptr: ptr; zbufptr, cuezbufptr, countbufptr: IntPtr; sumbufptr: LongPtr; str: string);
		var
			i, j, k,                                                 {loop indices}
			thispixel: integer;                              {the current pixel to be projected}
			offset, offsetinit: longint;                   {precomputed offsets into an image buffer}
			projaddr,                                            {buffer address for final projected image}
			opaaddr,                                              {buffer address for opacity surface projection component}
			brightcueaddr,                                   {buffer address for brightest-point projection with interior depth cues}
			zbufaddr,                                            {z-buffer address for surface projections (nearest-point/opacity)}
			cuezbufaddr,                                       {z-buffer address for brightest-point projection with interior depth cues}
			countbufaddr,                                     {buffer addr for counting pixels in mean-value projection}
			sumbufaddr: longint;                         {buffer addr for summing pixels in mean-value projection}
			z,                                                        {z-coordinate of points in current slice before rotation}
			ynew, znew,                                       {y- and z-coordinates of current point after rotation}
			zmax, zmin,                                       {z-coordinates of first and last slices before rotation}
			zmaxminuszmintimes100: longint;  {precomputed values to save time in loops}
			c100minusDepthCueInt, c100minusDepthCueSurf: integer;
			DepthCueIntlessthan100, DepthCueSurflessthan100: boolean;
			OpacityOrNearestPt, OpacityAndNotNearestPt: boolean;
			MeanVal, BrightestPt: boolean;
			ysintheta, ycostheta,                          {values used in rotational calculations before projection}
			zsintheta, zcostheta, ysinthetainit, ycosthetainit: longint;
			p,                                                        {pointer to final projected image buffer}
			op,                                                      {pointer to opacity surface projection buffer}
			bp: ptr;                                              {pointer to brightest-point projection buffer with interior depth cues}
			zp,                                                      {pointer to surface projection (nearest-point/opacity) z-buffer}
			qp,                                                      {pointer to z-buffer for brightest-point projection with interior depth cues}
			cp: IntPtr;                                          {pointer to buffer for counting pixels for mean-value projection}
			sp: LongPtr;                                       {pointer to mean-value summing buffer}
			width: integer;
			theLine: LineType;
	begin
		with BoundRect do begin                    {precompute values to avoid computations in loop}
				width := right - left;
				zmax := zcenter + projheight div 2; {find z-coordinates of first and last slices}
				zmin := zcenter - projheight div 2;
				zmaxminuszmintimes100 := 100 * (zmax - zmin);
				c100minusDepthCueInt := 100 - DepthCueInt;
				c100minusDepthCueSurf := 100 - DepthCueSurf;
				DepthCueIntlessthan100 := DepthCueInt < 100;
				DepthCueSurflessthan100 := DepthCueSurf < 100;
				OpacityOrNearestPt := (ProjectionMethod = NearestPoint) or (Opacity > 0);
				OpacityAndNotNearestPt := (Opacity > 0) and (ProjectionMethod <> NearestPoint);
				MeanVal := ProjectionMethod = MeanValue;
				BrightestPt := ProjectionMethod = BrightestPoint;
				projaddr := ord4(projptr);
				opaaddr := ord4(opaptr);
				brightcueaddr := ord4(brightcueptr);
				zbufaddr := ord4(zbufptr);
				cuezbufaddr := ord4(cuezbufptr);
				countbufaddr := ord4(countbufptr);
				sumbufaddr := ord4(sumbufptr);
				ycosthetainit := (top - ycenter - 1) * costheta;
				ysinthetainit := (top - ycenter - 1) * sintheta;
				offsetinit := ((projheight - bottom + top) div 2) * projwidth + (projwidth - right + left) div 2 - 1;
				for k := 1 to nSlices do begin
						SelectSlice(k);
						z := round((k - 1) * SliceInterval) - zcenter;
						zcostheta := z * costheta;
						zsintheta := z * sintheta;
						ycostheta := ycosthetainit;
						ysintheta := ysinthetainit;
						for j := top to bottom - 1 do begin               {read each row in the current slice}
								ycostheta := ycostheta + costheta;               {rotate about x-axis and find new y,z}
								ysintheta := ysintheta + sintheta;               {x-coordinates will not change}
								ynew := (ycostheta - zsintheta) div bigpowerof2 + ycenter - top;
								znew := (ysintheta + zcostheta) div bigpowerof2 + zcenter;
								offset := offsetinit + ynew * projwidth;
								GetLine(left, j, width, theLine);
								for i := 0 to width - 1 do begin             {read each pixel in current row and project it}
										thispixel := theLine[i];
										offset := offset + 1;
										if (offset >= ProjSize) or (offset < 0) then
											offset := 0;
										if (thispixel <= Upper) and (thispixel >= Lower) then begin
												if OpacityOrNearestPt then begin
														zp := IntPtr(zbufaddr + offset + offset);
														if (znew < zp^) then begin
																zp^ := znew;
																if OpacityAndNotNearestPt then begin
																		op := ptr(opaaddr + offset);
																		if (DepthCueSurflessthan100) then
																			op^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100)
																		else
																			op^ := thispixel;
																	end
																else begin
																		p := ptr(projaddr + offset);
																		if DepthCueSurflessthan100 then
																			p^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100)
																		else
																			p^ := thispixel;
																	end;
															end;
													end;     {NearestPoint case}
												if MeanVal then begin
														sp := LongPtr(sumbufaddr + offset + offset + offset + offset);
														sp^ := sp^ + thispixel;
														cp := IntPtr(countbufaddr + offset + offset);
														cp^ := cp^ + 1;
													end     {MeanValue case}
												else if BrightestPt then begin
														if (DepthCueIntlessthan100) then begin
																p := ptr(projaddr + offset);
																bp := ptr(brightcueaddr + offset);
																qp := IntPtr(cuezbufaddr + offset + offset);
																if (thispixel < BAND(bp^, 255)) or ((thispixel = BAND(bp^, 255)) and (znew < qp^)) then begin
																		bp^ := thispixel;                     {use z-buffer to ensure that if depth-cueing is on,  }
																		qp^ := znew;                            {the closer of two equally-bright points is displayed}
																		p^ := 255 - (DepthCueInt * (255 - thispixel) div 100 + (c100minusDepthCueInt) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100);
																	end;     {then}
															end
														else begin
																p := ptr(projaddr + offset);
																if (thispixel < BAND(p^, 255)) then
																	p^ := thispixel;
															end;     {else}
													end;     {BrightestPoint case}
											end;
									end;     {for j}
							end;     {for i}
						UpdateMeter(10 + (k * 90) div nSlices, str);
						if CommandPeriod then begin
								cancelled := true;
								beep;
								leave;
							end;
					end;     {for k}
			end;     {with}
	end;     {procedure DoOneProjectionX}



{******************************************************************************}
{*     This procedure projects each pixel of a volume (stack of slices) onto a plane as the volume rotates about    *}
{*  the y-axis.  Integer arithmetic, precomputation of values, and iterative addition rather than multiplication  *}
{*  inside a loop are used extensively to make the code run efficiently.  Projection parameters stored in global   *}
{*  variables determine how the projection will be performed.  This procedure returns various buffers which   *}
{*  are actually used by DoProject to find the final projected image for the volume of slices at the current angle.*}
{******************************************************************************}
	procedure DoOneProjectionY (nSlices, xcenter, zcenter: integer; projwidth, projheight: LongInt; costheta, sintheta: longint; projptr, opaptr, brightcueptr: ptr; zbufptr, cuezbufptr, countbufptr: IntPtr; sumbufptr: LongPtr; str: string);
		var
			i, j, k, thispixel: integer;
			offset, offsetinit: longint;
			projaddr, opaaddr, brightcueaddr, zbufaddr, cuezbufaddr, countbufaddr, sumbufaddr: longint;
			z, xnew, znew, zmax, zmin, zmaxminuszmintimes100: longint;
			c100minusDepthCueInt, c100minusDepthCueSurf: integer;
			DepthCueIntlessthan100, DepthCueSurflessthan100: boolean;
			OpacityOrNearestPt, OpacityAndNotNearestPt: boolean;
			MeanVal, BrightestPt: boolean;
			xsintheta, xcostheta, zsintheta, zcostheta, xsinthetainit, xcosthetainit: longint;
			p, op, bp: ptr;
			zp, qp, cp: IntPtr;
			sp: LongPtr;
			width: integer;
			aLine: LineType;
	begin
		with BoundRect do begin
				width := right - left;
				zmax := zcenter + projwidth div 2;
				zmin := zcenter - projwidth div 2;
				zmaxminuszmintimes100 := 100 * (zmax - zmin);
				c100minusDepthCueInt := 100 - DepthCueInt;
				c100minusDepthCueSurf := 100 - DepthCueSurf;
				DepthCueIntlessthan100 := DepthCueInt < 100;
				DepthCueSurflessthan100 := DepthCueSurf < 100;
				OpacityOrNearestPt := (ProjectionMethod = NearestPoint) or (Opacity > 0);
				OpacityAndNotNearestPt := (Opacity > 0) and (ProjectionMethod <> NearestPoint);
				MeanVal := ProjectionMethod = MeanValue;
				BrightestPt := ProjectionMethod = BrightestPoint;
				projaddr := ord4(projptr);
				opaaddr := ord4(opaptr);
				brightcueaddr := ord4(brightcueptr);
				zbufaddr := ord4(zbufptr);
				cuezbufaddr := ord4(cuezbufptr);
				countbufaddr := ord4(countbufptr);
				sumbufaddr := ord4(sumbufptr);
				xcosthetainit := (left - xcenter - 1) * costheta;
				xsinthetainit := (left - xcenter - 1) * sintheta;
				for k := 1 to nSlices do begin
						SelectSlice(k);
						z := round((k - 1) * SliceInterval) - zcenter;
						zcostheta := z * costheta;
						zsintheta := z * sintheta;
						offsetinit := ((projheight - bottom + top) div 2) * projwidth + (projwidth - right + left) div 2 - projwidth;
						for j := top to bottom - 1 do begin
								xcostheta := xcosthetainit;
								xsintheta := xsinthetainit;
								offsetinit := offsetinit + projwidth;
								GetLine(left, j, width, aLine);
								for i := 0 to width - 1 do begin
										thispixel := aLine[i];
										xcostheta := xcostheta + costheta;
										xsintheta := xsintheta + sintheta;
										if (thispixel <= Upper) and (thispixel >= Lower) then begin
												xnew := (xcostheta + zsintheta) div bigpowerof2 + xcenter - left;
												znew := (zcostheta - xsintheta) div bigpowerof2 + zcenter;
												offset := offsetinit + xnew;
												if (offset >= ProjSize) or (offset < 0) then
													offset := 0;
												if OpacityOrNearestPt then begin
														zp := IntPtr(zbufaddr + offset + offset);
														if (znew < zp^) then begin
																zp^ := znew;
																if OpacityAndNotNearestPt then begin
																		op := ptr(opaaddr + offset);
																		if (DepthCueSurflessthan100) then
																			op^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100)
																		else
																			op^ := thispixel;
																	end
																else begin
																		p := ptr(projaddr + offset);
																		if DepthCueSurflessthan100 then
																			p^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100)
																		else
																			p^ := thispixel;
																	end;
															end;
													end;     {NearestPoint case}
												if MeanVal then begin
														sp := LongPtr(sumbufaddr + offset + offset + offset + offset);
														sp^ := sp^ + thispixel;
														cp := IntPtr(countbufaddr + offset + offset);
														cp^ := cp^ + 1;
													end     {MeanValue case}
												else if BrightestPt then begin
														if (DepthCueIntlessthan100) then begin
																p := ptr(projaddr + offset);
																bp := ptr(brightcueaddr + offset);
																qp := IntPtr(cuezbufaddr + offset + offset);
																if (thispixel < BAND(bp^, 255)) or ((thispixel = BAND(bp^, 255)) and (znew < qp^)) then begin
																		bp^ := thispixel;
																		qp^ := znew;
																		p^ := 255 - (DepthCueInt * (255 - thispixel) div 100 + (c100minusDepthCueInt) * (255 - thispixel) * (zmax - znew) div zmaxminuszmintimes100);
																	end;     {then}
															end
														else begin
																p := ptr(projaddr + offset);
																if (thispixel < BAND(p^, 255)) then
																	p^ := thispixel;
															end;     {else}
													end;     {BrightestPoint case}
											end;
									end;     {for j}
							end;     {for i}
						UpdateMeter(10 + (k * 90) div nSlices, str);
						if CommandPeriod then begin
								cancelled := true;
								beep;
								leave;
							end;
					end;     {for k}
			end;     {with}
	end;     {procedure DoOneProjectionY}



{******************************************************************************}
{*     This procedure projects each pixel of a volume (stack of slices) onto a plane as the volume rotates about    *}
{*  the z-axis.  Integer arithmetic, precomputation of values, and iterative addition rather than multiplication  *}
{*  inside a loop are used extensively to make the code run efficiently.  Projection parameters stored in global   *}
{*  variables determine how the projection will be performed.  This procedure returns various buffers which   *}
{*  are actually used by DoProject to find the final projected image for the volume of slices at the current angle.*}
{******************************************************************************}
	procedure DoOneProjectionZ (nSlices, xcenter, ycenter, zcenter: integer; projwidth, projheight: LongInt; costheta, sintheta: longint; projptr, opaptr, brightcueptr: ptr; zbufptr, cuezbufptr, countbufptr: IntPtr; sumbufptr: LongPtr; str: string);
		var
			i, j, k, thispixel: integer;
			offset, offsetinit: longint;
			projaddr, opaaddr, brightcueaddr, zbufaddr, cuezbufaddr, countbufaddr, sumbufaddr: longint;
			z, xnew, ynew, zmax, zmin, zmaxminuszmintimes100: longint;
			c100minusDepthCueInt, c100minusDepthCueSurf: integer;
			DepthCueIntlessthan100, DepthCueSurflessthan100: boolean;
			OpacityOrNearestPt, OpacityAndNotNearestPt: boolean;
			MeanVal, BrightestPt: boolean;
			xsintheta, xcostheta, ysintheta, ycostheta: longint;
			xsinthetainit, xcosthetainit, ysinthetainit, ycosthetainit: longint;
			p, op, bp: ptr;
			zp, qp, cp: IntPtr;
			sp: LongPtr;
			width: integer;
			theLine: LineType;
	begin
		with BoundRect do begin
				width := right - left;
				zmax := zcenter + projwidth div 2;
				zmin := zcenter - projwidth div 2;
				zmaxminuszmintimes100 := 100 * (zmax - zmin);
				c100minusDepthCueInt := 100 - DepthCueInt;
				c100minusDepthCueSurf := 100 - DepthCueSurf;
				DepthCueIntlessthan100 := DepthCueInt < 100;
				DepthCueSurflessthan100 := DepthCueSurf < 100;
				OpacityOrNearestPt := (ProjectionMethod = NearestPoint) or (Opacity > 0);
				OpacityAndNotNearestPt := (Opacity > 0) and (ProjectionMethod <> NearestPoint);
				MeanVal := ProjectionMethod = MeanValue;
				BrightestPt := ProjectionMethod = BrightestPoint;
				projaddr := ord4(projptr);
				opaaddr := ord4(opaptr);
				brightcueaddr := ord4(brightcueptr);
				zbufaddr := ord4(zbufptr);
				cuezbufaddr := ord4(cuezbufptr);
				countbufaddr := ord4(countbufptr);
				sumbufaddr := ord4(sumbufptr);
				xcosthetainit := (left - xcenter - 1) * costheta;
				xsinthetainit := (left - xcenter - 1) * sintheta;
				ycosthetainit := (top - ycenter - 1) * costheta;
				ysinthetainit := (top - ycenter - 1) * sintheta;
				offsetinit := ((projheight - bottom + top) div 2) * projwidth + (projwidth - right + left) div 2 + left - 1;
				for k := 1 to nSlices do begin
						SelectSlice(k);
						z := round((k - 1) * SliceInterval) - zcenter;
						ycostheta := ycosthetainit;
						ysintheta := ysinthetainit;
						for j := top to bottom - 1 do begin
								ycostheta := ycostheta + costheta;
								ysintheta := ysintheta + sintheta;
								xcostheta := xcosthetainit;
								xsintheta := xsinthetainit;
								GetLine(left, j, width, theLine);
								for i := 0 to width - 1 do begin
										thisPixel := theLine[i];
										xcostheta := xcostheta + costheta;
										xsintheta := xsintheta + sintheta;
										if (thispixel <= Upper) and (thispixel >= Lower) then begin
												xnew := (xcostheta - ysintheta) div bigpowerof2 + xcenter - left;
												ynew := (xsintheta + ycostheta) div bigpowerof2 + ycenter - top;
												offset := offsetinit + ynew * projwidth + xnew;
												if (offset >= ProjSize) or (offset < 0) then
													offset := 0;
												if OpacityOrNearestPt then begin
														zp := IntPtr(zbufaddr + offset + offset);
														if (z < zp^) then begin
																zp^ := z;
																if OpacityAndNotNearestPt then begin
																		op := ptr(opaaddr + offset);
																		if (DepthCueSurflessthan100) then
																			op^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - z) div zmaxminuszmintimes100)
																		else
																			op^ := thispixel;
																	end
																else begin
																		p := ptr(projaddr + offset);
																		if DepthCueSurflessthan100 then
																			p^ := 255 - (DepthCueSurf * (255 - thispixel) div 100 + (c100minusDepthCueSurf) * (255 - thispixel) * (zmax - z) div zmaxminuszmintimes100)
																		else
																			p^ := thispixel;
																	end;
															end;
													end;     {NearestPoint case}
												if MeanVal then begin
														sp := LongPtr(sumbufaddr + offset + offset + offset + offset);
														sp^ := sp^ + thispixel;
														cp := IntPtr(countbufaddr + offset + offset);
														cp^ := cp^ + 1;
													end     {MeanValue case}
												else if BrightestPt then begin
														if (DepthCueIntlessthan100) then begin
																p := ptr(projaddr + offset);
																bp := ptr(brightcueaddr + offset);
																qp := IntPtr(cuezbufaddr + offset + offset);
																if (thispixel < BAND(bp^, 255)) or ((thispixel = BAND(bp^, 255)) and (z < qp^)) then begin
																		bp^ := thispixel;
																		qp^ := z;
																		p^ := 255 - (DepthCueInt * (255 - thispixel) div 100 + (c100minusDepthCueInt) * (255 - thispixel) * (zmax - z) div zmaxminuszmintimes100);
																	end;     {then}
															end
														else begin
																p := ptr(projaddr + offset);
																if (thispixel < BAND(p^, 255)) then
																	p^ := thispixel;
															end;     {else}
													end;     {BrightestPoint case}
											end;
									end;     {for j}
							end;     {for i}
						UpdateMeter(10 + (k * 90) div nSlices, str);
						if CommandPeriod then begin
								cancelled := true;
								beep;
								leave;
							end;
					end;     {for k}
			end;     {with}
	end;     {procedure DoOneProjectionZ}



{******************************************************************************}
{*     This code initializes buffers by filling them with uniform values.  The buffers may be filled with               *}
{*  one-byte, two-byte, or four-byte values (maybe others, too!).  Since multiple, huge buffers are                  *}
{*  allocated for projection calculations, we have decided to write this procedure in assembly language.               *}
{*  Assembly code written by Janice Keller.                                                                                                              *}
{******************************************************************************}
	procedure InitializeBuffer (p: ptr; size: longint; value, step: integer);
	inline
		$4E56, $0000,					{link  A6,#0}
		$48E7, $F880,          {movem.l  a0/d0-d4, -(sp)}
		$342E, $0004,         {move.w   4(a6),d2               step}
		$7000,                     {clr.l  d0                               set high word}
		$302E, $0006,         {move.w   6(a6),d0               value to set}
		$222E, $0008,         {move.l   8(a6),d1               projsize}
		$206E, $000C,         {movea.l  12(a6),a0             start address}
		$7600,                     {Test1 clr.l   d3                   set remainder to 0}
		$0881, $0000,        {bclr.l       #0,d1                  test for 1 and zero}
		$6702,                    {beq.b       Test2                    if 1, save it, else continue}
		$7601,                     {moveq.l     #1,d3                remainder is 1}
		$0881, $0001,        {Test2  bclr.l  #1,d1             test for 2 or 3}
		$6702,                    {beq.b  SetValues                   if found 1, set remainder}
		$5403,                    {addq.b     #2,d3                   remainder + 2}
		$7801,                    {SetValues moveq.l    #1,d4  decrement projsize by 1 for step 4}
		$0C02, $0004,        {cmp.b   #4,d2     if step is 4...}
		$6716,                    {beq.b     SetInit                    start initting}
		$7802,                    {moveq.l  #2,d4                   decr projsize by 2 for step 2}
		$0C02, $0002,        {cmp.b     #2,d2                   if step is 2...}
		$6708,                   {beq.b   DoubleIt                    just have to double}
		$7804,                    {moveq.l   #4,d4                  decre projsize by 4 for step 1}
		$1400,                    {move.b   d0,d2                    else we have a 1-er}
		$E148,                    {lsl.w   #8,d0                       lo to hi}
		$1002,                    {move.b  d2,d0                     reset lo}
		$3400,                   {DoubleIt  move.w d0,d2        save the lo word}
		$4840,                   {swap  d0                              move lo to hi}
		$3002,                   {move.w  d2,d0                     reset lo}
		$20C0,                   {SetInit   move.l   d0,(a0)+    set this address}
		$9284,                   {sub.l  d4,d1                         decr projsize}
		$0C81, $0000, $0000,  {cmp.l   #0,d1               are we done yet}
		$6EF4,                   {bgt.b  SetInit												no}
		$0C03, $0000,      {Remainder  cmp.b  #0,d3     is there anything left}
		$672C,                  {beq.b   Exit                           no, all done}
		$5383,                  {subq.l  #1,d3			                  0-2, not 1-3}
		$302E, $0006,      {move.w  6(a6),d0                 get the value}
		$342E, $0004,      {move.w 4(a6),d2                  get the step}
		$0C02, $0004,        {cmp.b  #4,d2                        is this a step 4}
		$6608,                     {bne.b   Teststep2                  no}
		$20C0,                   {Loop4   move.l  d0,(a0)+      yes, set long}
		$51CB, $FFFC,       {dbra   d3,Loop4                    next one}
		$6014,                  {bra.b  Exit}
		$0C02, $0002,      {TestStep2   cmp.b  #2,d2     is this a step 2}
		$6608,                 {bne.b   Loop1                       no, must be a 1}
		$30C0,                 {Loop2  move.w  d0,(a0)+     yes, set word}
		$51CB, $FFFC,       {dbra   d3,Loop2                    next one}
		$6006,                 {bra.b  Exit}
		$10C0,                 {Loop1  move.b   d0,(a0)+      set bytes}
		$51CB, $FFFC,      {dbra  d3,Loop1                     next one}
		$4CDF, $011F,       {Exit  movem.l    (sp)+,a0/d0-d4}
		$4E5E,                   {unlk    a6}
		$DEFC, $000C;      {add.w   #12,sp}
{this Pascal code works except for the fact that the pointer must change type with step}
{so if you want to use the Pascal, you'll need a case statement for each allowable size}
{$IFC false}
	var
		offset: longint;
begin
	for offset := 0 to size - 1 do begin
			p^ := value;
			p := Ptr(ord4(p) + step);
		end;     {for}
end;     {procedure InitializeBuffer}
{$ENDC}


{******************************************************************************}
{*     This procedure creates a sequence of projections of a rotating volume (stack of slices) onto a plane using   *}
{*  nearest-point (surface), brightest-point, or mean-value projection or a weighted combination of nearest- *}
{*  point projection with either of the other two methods (partial opacity).  The user may choose to rotate the   *}
{*  volume about any of the three orthogonal axes (x,y, or z), make portions of the volume transparent, or add  *}
{*  a greater degree of visual realism by employing depth cues.                                                                               *}
{******************************************************************************}
procedure DoProjection;
	var
		tport: Grafptr;
		nSlices: integer;                                      {number of slices in volume}
		projwidth, projheight: LongInt;              {dimensions of projection image}
		xcenter, ycenter, zcenter: integer;         {coordinates of center of volume of rotation}
		theta: integer;                                          {current angle of rotation in degrees}
		thetarad: real;                                          {current angle of rotation in radians}
		sintheta, costheta: longint;                      {sine and cosine of current angle}
		xsinthetainit, xcosthetainit: longint;      {precomputed products to avoid mult in loops}
		offset, MemoryRequired: longint;
		p, op, bp, projptr, opaptr, brightcueptr: ptr;
		zp, zbufptr, qp, cuezbufptr, countbufptr, cp: IntPtr;
		sp, sumbufptr: LongPtr;
		curval, prevval, nextval, aboveval, belowval: integer;
		ignore: integer;                                        {irrelevant return value from a function}
		ticksinit, tickstogo, TicksForOneProjection: longint;
		str, TimeStr, seconds: str255;
		SaveProjectionsTemp, AutoSelectAll, AllocatingBuffers: boolean;
		n, nProjections, angle: integer;
		SourceInfo, DestInfo: InfoPtr;

	procedure Abort;
	begin
		DisposeProjectionPtrs(projptr, opaptr, brightcueptr, zbufptr, countbufptr, cuezbufptr, sumbufptr);
		if AllocatingBuffers and (MaxBlock > 20000) then
			PutMessage('Insufficient Memory.');
		macro := false;
		exit(DoProjection);
	end;

begin
	ShowWatch;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	if NotInBounds then
		exit(DoProjection);
	cancelled := false;
	SourceInfo := Info;
	GetPort(tPort);
	with Info^ do begin
			SetPort(GrafPtr(osPort));
			BoundRect := Roirect;
		end;
	if (AngleInc = 0) and (TotalAngle <> 0) then
		AngleInc := 5;
	angle := 0;
	nProjections := 0;
	if AngleInc = 0 then
		nProjections := 1
	else
		while angle <= TotalAngle do begin
				nProjections := nProjections + 1;
				angle := angle + AngleInc;
			end;
	if angle > 360 then
		nProjections := nProjections - 1;
	if nProjections <= 0 then
		nProjections := 1;
	nSlices := Info^.StackInfo^.nSlices;         {get number of slices in volume}
	with BoundRect do begin
			xcenter := (left + right) div 2;          {find center of volume of rotation}
			ycenter := (top + bottom) div 2;
			zcenter := round(nSlices * SliceInterval / 2.0);
			if MinProjSize and (AxisOfRotation <> ZAxis) then begin
					case AxisOfRotation of                    {find dimensions of projection image}
						XAxis:  begin
								projheight := round(sqrt(sqr(nSlices * SliceInterval) + longint(bottom - top) * (bottom - top)));
								projwidth := right - left;
							end;     {XAxis}
						YAxis:  begin
								projwidth := round(sqrt(sqr(nSlices * SliceInterval) + longint(right - left) * (right - left)));
								projheight := bottom - top;
							end;     {YAxis}
					end;     {case}
				end     {then}
			else begin
					projwidth := round(sqrt(sqr(nSlices * SliceInterval) + longint(right - left) * (right - left)));
					projheight := round(sqrt(sqr(nSlices * SliceInterval) + longint(bottom - top) * (bottom - top)));
				end;     {else make all windows the same size regardless of rotation axis}
		end;     {with BoundRect}
	if odd(projwidth) then
		projwidth := projwidth + 1;
	projptr := nil;
	zbufptr := nil;
	opaptr := nil;
	brightcueptr := nil;
	cuezbufptr := nil;
	sumbufptr := nil;
	countbufptr := nil;
	AllocatingBuffers := true;
	projsize := projwidth * projheight;
	projptr := NewPtr(projsize);
	if projptr = nil then
		Abort;
	if (ProjectionMethod = NearestPoint) or (Opacity > 0) then begin    {get other required buffers}
			zbufptr := IntPtr(NewPtr(projsize * 2));
			if zbufptr = nil then
				Abort;
		end;
	if (Opacity > 0) and (ProjectionMethod <> NearestPoint) then begin
			opaptr := NewPtr(projsize);
			if opaptr = nil then
				Abort;
		end;
	if (ProjectionMethod = BrightestPoint) and (DepthCueInt < 100) then begin
			brightcueptr := NewPtr(projsize);
			cuezbufptr := IntPtr(NewPtr(projsize * 2));
			if (brightcueptr = nil) or (cuezbufptr = nil) then
				abort;
		end;
	if (ProjectionMethod = MeanValue) then begin
			sumbufptr := LongPtr(NewPtr(projsize * 4));
			countbufptr := IntPtr(NewPtr(projsize * 2));
			if (sumbufptr = nil) or (countbufptr = nil) then
				abort;
		end;
	AllocatingBuffers := false;
	SaveProjectionsTemp := FALSE;              {check whether we have enough memory to open}
	MemoryRequired := nProjections * projsize + SizeOf(PicInfo) + MinFree;
	if (MemoryRequired > FreeMem) and not (SaveProjections) then begin
			str := 'Insufficient memory to create entire stack of projections.  Projection images will be saved to disk.';
			if (PutMessageWithCancel(str) = cancel) then
				Abort;
			SaveProjections := TRUE;
			SaveProjectionsTemp := TRUE;
		end;
	if (SaveProjections) then begin           {prepare to save projections as created if desired}
			SaveAsWhat := AsTiff;
			SaveAllState := SaveAllStage1;
		end;
	TimeStr := '?';
	theta := InitAngle;                                     {begin rotation with user-specified angle}
	ticksinit := TickCount;
	for n := 1 to nProjections do begin
			if (SaveProjections) or (n = 1) then begin    {open new window or stack slice}
					if SaveProjections then
						case AxisOfRotation of
							XAxis: 
								str := concat('Projection X ', long2str(theta));
							YAxis: 
								str := concat('Projection Y ', long2str(theta));
							ZAxis: 
								str := concat('Projection Z ', long2str(theta));
						end
					else
						str := 'Projections';
					if not NewPicWindow(str, projwidth, projheight) then
						Abort;
				end
			else if not (AddSlice(false)) then
				Abort;
			str := concat('Projecting: ', long2str(n), '/', long2str(nProjections), ' (', long2str(Theta), '¡', ', ', TimeStr, ')');
			ShowMeter;
			UpdateMeter(0, str);
			thetarad := theta * pi / 180;
			costheta := round(bigpowerof2 * cos(thetarad));
			sintheta := round(bigpowerof2 * sin(thetarad));
			p := projptr;                                          {initialize required projection buffers}
			InitializeBuffer(p, projsize, 255, 1);
			if (ProjectionMethod = NearestPoint) or (Opacity > 0) then begin
					zp := zbufptr;
					InitializeBuffer(Ptr(zp), projsize, 32767, 2);
				end;     {then}
			if (Opacity > 0) and (ProjectionMethod <> NearestPoint) then begin
					op := opaptr;
					InitializeBuffer(op, projsize, 255, 1);
				end;     {then}
			if (ProjectionMethod = MeanValue) then begin
					sp := sumbufptr;
					cp := countbufptr;
					InitializeBuffer(Ptr(sp), projsize, 0, 4);
					InitializeBuffer(Ptr(cp), projsize, 0, 2);
				end;
			if (ProjectionMethod = BrightestPoint) and (DepthCueInt < 100) then begin
					bp := brightcueptr;
					InitializeBuffer(bp, projsize, 255, 1);
					qp := cuezbufptr;
					InitializeBuffer(Ptr(qp), projsize, 255, 2);
				end;     {then}
			UpdateMeter(10, str);
			DestInfo := Info;
			Info := SourceInfo;
			case AxisOfRotation of
				XAxis: 
					DoOneProjectionX(nSlices, ycenter, zcenter, projwidth, projheight, costheta, sintheta, projptr, opaptr, brightcueptr, zbufptr, cuezbufptr, countbufptr, sumbufptr, str);
				YAxis: 
					DoOneProjectionY(nSlices, xcenter, zcenter, projwidth, projheight, costheta, sintheta, projptr, opaptr, brightcueptr, zbufptr, cuezbufptr, countbufptr, sumbufptr, str);
				ZAxis: 
					DoOneProjectionZ(nSlices, xcenter, ycenter, zcenter, projwidth, projheight, costheta, sintheta, projptr, opaptr, brightcueptr, zbufptr, cuezbufptr, countbufptr, sumbufptr, str);
			end;
			Info := DestInfo;
			if n = 1 then
				TicksForOneProjection := TickCount - TicksInit;
			TicksToGo := (nProjections - n) * TicksForOneProjection;
			NumToString((TicksToGo div 60) mod 60, seconds);
			if length(seconds) = 1 then
				seconds := concat('0', seconds);
			timestr := concat(long2str(TicksToGo div 3600), ':', seconds);
			if (ProjectionMethod = MeanValue) then begin
					p := projptr;                                    {calculate the mean-value image from returned info}
					sp := sumbufptr;
					cp := countbufptr;
					for offset := 0 to projsize - 1 do begin
							if (cp^ > 0) then
								p^ := sp^ div cp^;
							p := ptr(ord4(p) + 1);
							sp := LongPtr(ord4(sp) + 4);
							cp := IntPtr(ord4(cp) + 2);
						end;     {for}
				end;     {then}
			if (Opacity > 0) and (ProjectionMethod <> NearestPoint) then begin
					p := projptr;                                    {calculate surface proj component (opacity on)}
					op := opaptr;                                     {and combine with another proj component}
					for offset := 0 to projsize - 1 do begin
							p^ := (Opacity * BAND(op^, 255) + (100 - Opacity) * BAND(p^, 255)) div 100;
							p := ptr(ord4(p) + 1);
							op := ptr(ord4(op) + 1);
						end;     {for}
				end;     {then}
			if (AxisOfRotation = ZAxis) then begin  {interpolate for z-axis rotation}
					p := ptr(ord4(projptr) + projwidth);
					for offset := projwidth to projsize - 1 - projwidth do begin
							curval := BAND(p^, 255);
							prevval := BAND(ptr(ord4(p) - 1)^, 255);
							nextval := BAND(ptr(ord4(p) + 1)^, 255);
							aboveval := BAND(ptr(ord4(p) - projwidth)^, 255);
							belowval := BAND(ptr(ord4(p) + projwidth)^, 255);
							if (curval = 255) and (prevval <> 255) and (nextval <> 255) and (aboveval <> 255) and (belowval <> 255) then
								p^ := (prevval + nextval + aboveval + belowval) div 4;
							p := ptr(ord4(p) + 1);
						end;
				end;
			if (SaveProjections) or (n = 1) then
				BlockMove(projptr, Info^.PicBaseAddr, projsize)         {whole ROI write to projection image}
			else
				BlockMove(projptr, Info^.StackInfo^.PicBaseH[n]^, projsize);
			UpdateMeter(-1, '');        {dispose of meter}
			if cancelled then begin
					if n > 1 then
						DeleteSlice;
					leave;
				end;
			if (SaveProjections) then begin
					SaveAs('', 0);                                    {just save and put away current image after creating it}
					ignore := CloseAWindow(info^.wptr);
				end
			else if n = 1 then begin  {create new stack for first projection if not saving projections}
					if not MakeStackFromWindow then
						Abort
				end;
			theta := (theta + AngleInc) mod 360;
			UpdatePicWindow;
		end;     {for}
	SaveAllState := NoSaveAll;
	if SaveProjectionsTemp then                 {turn this back off if we turned it on due to lack of memory}
		SaveProjections := FALSE;
	DisposeProjectionPtrs(projptr, opaptr, brightcueptr, zbufptr, countbufptr, cuezbufptr, sumbufptr);
	SetPort(tPort);
	DestInfo := info;
	info := SourceInfo;
	SelectSlice(info^.StackInfo^.CurrentSlice);
	if AutoSelectAll then
		KillRoi;
	info := DestInfo;
end;     {procedure DoProjection}


{******************************************************************************}
{*     This procedure allows the user to set parameters for projection in a large dialog box.                                  *}
{******************************************************************************}
function ShowProjectDialogBox: boolean;
	const
		ProjectOptionsID = 128;
		SliceIntervalID = 3;
		InitAngleID = 4;
		TotalAngleID = 5;
		AngleIncID = 6;
		TransparencyLowerID = 7;
		TransparencyUpperID = 8;
		OpacityID = 9;
		DepthCueSurfID = 10;
		DepthCueIntID = 11;
		RotationAboutXID = 12;
		RotationAboutYID = 13;
		RotationAboutZID = 14;
		SaveProjectionsID = 15;
		MinProjSizeID = 16;
		NearestID = 28;
		BrightestID = 29;
		MeanID = 30;
	var
		mylog: Dialogptr;                                           {pointer to dialog box}
		i, item, alldone: integer;
		SaveInitAngle, SaveTotalAngle, SaveAngleInc: integer;
		SaveOpacity: integer;
		SaveAxisOfRotation: AxisType;
		SaveSaveProjections, SaveCloseSlices, SaveMinProjSize: Boolean;
		PercentSurf, PercentInt: integer;
begin
	InitCursor;
	SliceInterval := info^.StackInfo^.SliceSpacing;
	if SliceInterval <= 0.0 then
		SliceInterval := 1.0;
	SaveInitAngle := InitAngle;
	SaveTotalAngle := TotalAngle;
	SaveAngleInc := AngleInc;
	SaveOpacity := Opacity;
	SaveAxisOfRotation := AxisOfRotation;
	SaveSaveProjections := SaveProjections;
	SaveMinProjSize := MinProjSize;
	PercentSurf := 100 - DepthCueSurf;
	PercentInt := 100 - DepthCueInt;
	if DensitySlicing then
		with info^ do begin
				lower := SliceStart;
				upper := SliceEnd;
			end
	else begin
			lower := TransparencyLower;
			upper := TransparencyUpper;
		end;
	mylog := GetNewDialog(ProjectOptionsID, nil, pointer(-1));
	SetDReal(MyLog, SliceIntervalID, SliceInterval, 1);
	SelIText(MyLog, SliceIntervalID, 0, 32767);
	SetDNum(MyLog, InitAngleID, InitAngle);
	SetDNum(MyLog, TotalAngleID, TotalAngle);
	SetDNum(MyLog, AngleIncID, AngleInc);
	SetDNum(MyLog, TransparencyLowerID, lower);
	SetDNum(MyLog, TransparencyUpperID, upper);
	SetDNum(MyLog, OpacityID, Opacity);
	SetDNum(MyLog, DepthCueSurfID, PercentSurf);
	SetDNum(MyLog, DepthCueIntID, PercentInt);
	OutlineButton(MyLog, ok, 16);
	SetDialogItem(MyLog, RotationAboutXID, ord(AxisOfRotation = XAxis));
	SetDialogItem(MyLog, RotationAboutYID, ord(AxisOfRotation = YAxis));
	SetDialogItem(MyLog, RotationAboutZID, ord(AxisOfRotation = ZAxis));
	SetDialogItem(MyLog, NearestID, ord(ProjectionMethod = NearestPoint));
	SetDialogItem(MyLog, BrightestID, ord(ProjectionMethod = BrightestPoint));
	SetDialogItem(MyLog, MeanID, ord(ProjectionMethod = MeanValue));
	SetDialogItem(MyLog, SaveProjectionsID, ord(SaveProjections));
	SetDialogItem(MyLog, MinProjSizeID, ord(MinProjSize));
	alldone := 0;
	repeat  {if we don't do it this way, ModalDialog throws us into code checking after each keystroke}
		repeat   {meaning you can't type in a 2 digit number}
			ModalDialog(nil, item);
			if item = SaveProjectionsID then begin
					SaveProjections := not SaveProjections;
					SetDialogItem(MyLog, SaveProjectionsID, ord(SaveProjections));
				end;
			if item = MinProjSizeID then begin
					MinProjSize := not MinProjSize;
					SetDialogItem(MyLog, MinProjSizeID, ord(MinProjSize));
				end;
			if (item = RotationAboutXID) or (item = RotationAboutYID) or (item = RotationAboutZID) then begin
					case item of
						RotationAboutXID: 
							AxisOfRotation := XAxis;
						RotationAboutYID: 
							AxisOfRotation := YAxis;
						RotationAboutZID: 
							AxisOfRotation := ZAxis;
					end;     {case}
					SetDialogItem(MyLog, RotationAboutXID, ord(AxisOfRotation = XAxis));
					SetDialogItem(MyLog, RotationAboutYID, ord(AxisOfRotation = YAxis));
					SetDialogItem(MyLog, RotationAboutZID, ord(AxisOfRotation = ZAxis));
				end;
			if (item >= nearestID) and (item <= MeanID) then begin
					case item of
						NearestID: 
							ProjectionMethod := NearestPoint;
						BrightestID: 
							ProjectionMethod := BrightestPoint;
						MeanID: 
							ProjectionMethod := MeanValue;
					end;
					SetDialogItem(MyLog, NearestID, ord(projectionMethod = NearestPoint));
					SetDialogItem(MyLog, BrightestID, ord(projectionMethod = BrightestPoint));
					SetDialogItem(MyLog, MeanID, ord(projectionMethod = MeanValue));
				end;
		until (item = ok) or (item = cancel);
		alldone := 1;
		if (item = ok) then begin  {check all the values that could have been entered, don't know which were changed}
				SliceInterval := GetDReal(MyLog, SliceIntervalID);
				if (SliceInterval <= 0.0) or (SliceInterval > 1023.0) then begin
						SliceInterval := info^.StackInfo^.SliceSpacing;
						SetDReal(MyLog, SliceIntervalID, SliceInterval, 1);
						beep;
						alldone := 0;
					end;  {if SliceInterval}
				InitAngle := GetDNum(MyLog, InitAngleID);
				if (InitAngle < 0) or (InitAngle > 360) then begin
						InitAngle := SaveInitAngle;
						SetDNum(MyLog, InitAngleID, InitAngle);
						beep;
						alldone := 0;
					end;  {if InitAngle}
				TotalAngle := GetDNum(MyLog, TotalAngleID);
				if (TotalAngle < 0) or (TotalAngle > 360) then begin
						TotalAngle := SaveTotalAngle;
						SetDNum(MyLog, TotalAngleID, TotalAngle);
						beep;
						alldone := 0;
					end;  {if TotalAngle}
				AngleInc := GetDNum(MyLog, AngleIncID);
				if (AngleInc < 0) or (AngleInc > 360) then begin
						AngleInc := SaveAngleInc;
						SetDNum(MyLog, AngleIncID, AngleInc);
						beep;
						alldone := 0;
					end;  {if AngleInc}
				lower := GetDNum(MyLog, TransparencyLowerID);
				if (lower < 0) or (lower > 255) then begin
						lower := TransparencyLower;
						SetDNum(MyLog, TransparencyLowerID, lower);
						beep;
						alldone := 0;
					end;  {if TransparencyLower}
				upper := GetDNum(MyLog, TransparencyUpperID);
				if (upper < 0) or (upper > 255) then begin
						upper := TransparencyUpper;
						SetDNum(MyLog, TransparencyUpperID, upper);
						beep;
						alldone := 0;
					end;  {if TransparencyUpper}
				Opacity := GetDNum(MyLog, OpacityID);
				if (Opacity < 0) or (Opacity > 100) then begin
						Opacity := SaveOpacity;
						SetDNum(MyLog, OpacityID, Opacity);
						beep;
						alldone := 0;
					end;  {if Opacity}
				PercentSurf := GetDNum(MyLog, DepthCueSurfID);
				if (PercentSurf < 0) or (PercentSurf > 100) then begin
						PercentSurf := 100 - DepthCueSurf;
						SetDNum(MyLog, DepthCueSurfID, PercentSurf);
						beep;
						alldone := 0;
					end;  {if DepthCueSurf}
				PercentInt := GetDNum(MyLog, DepthCueIntID);
				if (PercentInt < 0) or (PercentInt > 100) then begin
						PercentInt := 100 - DepthCueInt;
						SetDNum(MyLog, DepthCueIntID, PercentInt);
						beep;
						alldone := 0;
					end;  {if DepthCueInt}
				info^.StackInfo^.SliceSpacing := SliceInterval;
			end;
	until (alldone = 1);
	DisposDialog(mylog);
	if item = cancel then begin                       {if Cancel, keep the saved values}
			InitAngle := SaveInitAngle;
			TotalAngle := SaveTotalAngle;
			AngleInc := SaveAngleInc;
			Opacity := SaveOpacity;
			AxisOfRotation := SaveAxisOfRotation;
			SaveProjections := SaveSaveProjections;
			MinProjSize := SaveMinProjSize;
			ShowProjectDialogBox := false;
		end
	else begin
			if not DensitySlicing then begin
					TransparencyLower := lower;
					TransparencyUpper := upper;
				end;
			DepthCueSurf := 100 - PercentSurf;
			DepthCueInt := 100 - PercentInt;
			ShowProjectDialogBox := true;
		end;
end;


procedure Project;
begin
	if ShowProjectDialogBox then
		DoProjection;
end;


end.